"""
    ForwardOperators

Comprehensive forward operator suite for computing H(x) - the mapping from model state 
space to observation space for all atmospheric observation types used in GSI. This module
implements sophisticated interpolation algorithms, coordinate transformations, and 
specialized operators for different observation platforms.

# Mathematical Framework

The forward operator H maps the model state vector x to observation space:
```
y = H(x) + ε
```

Where:
- `x`: Model state vector (temperature, humidity, winds, etc.)
- `H(x)`: Forward operator mapping state to observation space
- `y`: Observations in their measurement space
- `ε`: Observation error

For linearized data assimilation, we also need:
- **Tangent Linear**: `δy = H'(x) δx` - linearized forward operator  
- **Adjoint**: `δx* = H'ᴵ(x) δy*` - transpose of linearized operator

# Observation Operator Types

## Conventional Observations
- **Temperature**: Spatial interpolation with pressure coordinate conversion
- **Humidity**: Interpolation with saturation checks and unit conversions
- **Winds**: Vector interpolation with proper coordinate transformation  
- **Pressure**: Hydrostatic adjustment and terrain correction

## Satellite Observations  
- **Radiances**: Radiative transfer model interface (CRTM)
- **Satellite Winds**: Height assignment and quality control
- **Retrievals**: Multi-channel regression and error propagation

## Specialized Observations
- **GPS Radio Occultation**: Ray path integration and Abel transform
- **Weather Radar**: Reflectivity/velocity forward modeling
- **Lightning**: Flash rate parameterization from convective indicators
- **Aerosols**: Optical property computation from mass concentrations

# Interpolation Methods

Multiple interpolation schemes optimized for different use cases:
- **Bilinear**: Fast interpolation for regular grids
- **Bicubic**: Higher accuracy for smooth fields  
- **Spline**: Maximum smoothness for sparse observations
- **Conservative**: Mass-conserving interpolation for integrated quantities
- **Distance-weighted**: Flexible interpolation for irregular grids
"""
module ForwardOperators

using LinearAlgebra
using Interpolations  
using StaticArrays
using ..ObservationTypes
using ..GSICoreAnalysis: AbstractAnalysisConfig, AbstractStateVector
using ..StateVectors: StateVector
using ..GridOperations: GridConfig, geographic_to_grid_coordinates

# Export main types and functions
export AbstractForwardOperator, LinearForwardOperator, NonlinearForwardOperator
export TemperatureOperator, HumidityOperator, WindOperator, PressureOperator
export RadianceOperator, GPSROOperator, RadarOperator
export InterpolationMethod, BilinearInterpolation, BicubicInterpolation, SplineInterpolation
export apply_forward_operator, apply_tangent_linear, apply_adjoint
export compute_observation_equivalent, interpolate_to_observation_location
export ForwardOperatorConfig, create_operator_suite

"""
    InterpolationMethod

Enumeration of available spatial interpolation methods.
"""
@enum InterpolationMethod begin
    BILINEAR_INTERPOLATION
    BICUBIC_INTERPOLATION
    SPLINE_INTERPOLATION
    CONSERVATIVE_INTERPOLATION
    DISTANCE_WEIGHTED_INTERPOLATION
    NEAREST_NEIGHBOR_INTERPOLATION
end

"""
    AbstractForwardOperator{T<:AbstractFloat}

Abstract base type for all forward operators.
"""
abstract type AbstractForwardOperator{T<:AbstractFloat} end

"""
    LinearForwardOperator{T<:AbstractFloat}

Base type for linear forward operators where H(x₁ + x₂) = H(x₁) + H(x₂).
"""
abstract type LinearForwardOperator{T<:AbstractFloat} <: AbstractForwardOperator{T} end

"""
    NonlinearForwardOperator{T<:AbstractFloat}

Base type for nonlinear forward operators requiring linearization.
"""
abstract type NonlinearForwardOperator{T<:AbstractFloat} <: AbstractForwardOperator{T} end

"""
    ForwardOperatorConfig{T<:AbstractFloat}

Configuration for forward operator computation.

# Fields
- `interpolation_method::InterpolationMethod`: Default interpolation method
- `boundary_treatment::Symbol`: How to handle domain boundaries (:extrapolate, :missing, :nearest)
- `vertical_coordinate::Symbol`: Vertical coordinate system (:pressure, :height, :sigma)
- `quality_control_checks::Bool`: Enable QC checks during operator application
- `parallel_processing::Bool`: Enable parallel processing for multiple observations
- `cache_interpolation_weights::Bool`: Cache interpolation weights for repeated use
- `tolerance::T`: Numerical tolerance for convergence checks
- `max_iterations::Int`: Maximum iterations for nonlinear operators
"""
struct ForwardOperatorConfig{T<:AbstractFloat}
    interpolation_method::InterpolationMethod
    boundary_treatment::Symbol
    vertical_coordinate::Symbol
    quality_control_checks::Bool
    parallel_processing::Bool
    cache_interpolation_weights::Bool
    tolerance::T
    max_iterations::Int
    
    function ForwardOperatorConfig{T}() where T
        new{T}(
            BILINEAR_INTERPOLATION,  # interpolation_method
            :extrapolate,            # boundary_treatment  
            :pressure,               # vertical_coordinate
            true,                    # quality_control_checks
            false,                   # parallel_processing
            true,                    # cache_interpolation_weights
            T(1e-6),                # tolerance
            50                       # max_iterations
        )
    end
end

ForwardOperatorConfig() = ForwardOperatorConfig{Float64}()

"""
    InterpolationWeights{T<:AbstractFloat}

Pre-computed interpolation weights for efficient repeated interpolation.

# Fields  
- `grid_indices::Vector{CartesianIndex{3}}`: Grid point indices for interpolation
- `weights::Vector{T}`: Interpolation weights corresponding to grid points
- `is_valid::Bool`: Whether the interpolation point is within the valid domain
"""
struct InterpolationWeights{T<:AbstractFloat}
    grid_indices::Vector{CartesianIndex{3}}
    weights::Vector{T}
    is_valid::Bool
end

"""
    TemperatureOperator{T<:AbstractFloat}

Forward operator for temperature observations.

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration
- `cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}`: Cached interpolation weights
- `pressure_levels::Vector{T}`: Model pressure levels
- `apply_lapse_rate_correction::Bool`: Apply lapse rate correction for elevation differences
"""
mutable struct TemperatureOperator{T<:AbstractFloat} <: LinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}
    pressure_levels::Vector{T}
    apply_lapse_rate_correction::Bool
    
    function TemperatureOperator{T}(config::ForwardOperatorConfig{T},
                                   pressure_levels::Vector{T}) where T
        new{T}(
            config,
            Dict{Tuple{T,T,T},InterpolationWeights{T}}(),
            pressure_levels,
            true
        )
    end
end

"""
    HumidityOperator{T<:AbstractFloat}

Forward operator for humidity observations (specific, relative, mixing ratio).

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration
- `cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}`: Cached interpolation weights
- `humidity_type::Symbol`: Type of humidity (:specific, :relative, :mixing_ratio)
- `saturation_check::Bool`: Enforce physical saturation limits
- `temperature_field::Union{Array{T,3},Nothing}`: Temperature field for humidity conversions
"""
mutable struct HumidityOperator{T<:AbstractFloat} <: NonlinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}
    humidity_type::Symbol
    saturation_check::Bool
    temperature_field::Union{Array{T,3},Nothing}
    
    function HumidityOperator{T}(config::ForwardOperatorConfig{T},
                                humidity_type::Symbol = :specific) where T
        new{T}(
            config,
            Dict{Tuple{T,T,T},InterpolationWeights{T}}(),
            humidity_type,
            true,
            nothing
        )
    end
end

"""
    WindOperator{T<:AbstractFloat}

Forward operator for wind component observations.

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration
- `cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}`: Cached interpolation weights
- `wind_component::Symbol`: Wind component (:u, :v, :speed, :direction)
- `coordinate_system::Symbol`: Coordinate system (:geographic, :grid, :natural)
- `rotation_angles::Union{Array{T,2},Nothing}`: Grid rotation angles for coordinate transformation
"""
mutable struct WindOperator{T<:AbstractFloat} <: LinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}
    wind_component::Symbol
    coordinate_system::Symbol
    rotation_angles::Union{Array{T,2},Nothing}
    
    function WindOperator{T}(config::ForwardOperatorConfig{T},
                            wind_component::Symbol = :u,
                            coordinate_system::Symbol = :geographic) where T
        new{T}(
            config,
            Dict{Tuple{T,T,T},InterpolationWeights{T}}(),
            wind_component,
            coordinate_system,
            nothing
        )
    end
end

"""
    PressureOperator{T<:AbstractFloat}

Forward operator for pressure observations.

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration
- `cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}`: Cached interpolation weights
- `pressure_type::Symbol`: Pressure type (:surface, :sea_level, :station)
- `terrain_elevation::Union{Array{T,2},Nothing}`: Model terrain elevation
- `reference_atmosphere::Union{Vector{T},Nothing}`: Reference atmosphere profile
"""
mutable struct PressureOperator{T<:AbstractFloat} <: LinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    cached_weights::Dict{Tuple{T,T,T},InterpolationWeights{T}}
    pressure_type::Symbol
    terrain_elevation::Union{Array{T,2},Nothing}
    reference_atmosphere::Union{Vector{T},Nothing}
    
    function PressureOperator{T}(config::ForwardOperatorConfig{T},
                                pressure_type::Symbol = :surface) where T
        new{T}(
            config,
            Dict{Tuple{T,T,T},InterpolationWeights{T}}(),
            pressure_type,
            nothing,
            nothing
        )
    end
end

"""
    RadianceOperator{T<:AbstractFloat}

Forward operator for satellite radiance observations using radiative transfer models.

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration
- `satellite_id::String`: Satellite platform identifier
- `instrument_id::String`: Instrument identifier  
- `channel_list::Vector{Int}`: List of channels to process
- `crtm_interface::Union{Any,Nothing}`: Interface to CRTM radiative transfer model
- `surface_type_map::Union{Array{Int,2},Nothing}`: Surface type classification
- `cloud_detection::Bool`: Enable cloud detection and handling
"""
mutable struct RadianceOperator{T<:AbstractFloat} <: NonlinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    satellite_id::String
    instrument_id::String
    channel_list::Vector{Int}
    crtm_interface::Union{Any,Nothing}
    surface_type_map::Union{Array{Int,2},Nothing}
    cloud_detection::Bool
    
    function RadianceOperator{T}(config::ForwardOperatorConfig{T},
                                satellite::String, instrument::String,
                                channels::Vector{Int}) where T
        new{T}(
            config, satellite, instrument, channels,
            nothing, nothing, true
        )
    end
end

"""
    GPSROOperator{T<:AbstractFloat}

Forward operator for GPS Radio Occultation observations.

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration
- `observation_type::Symbol`: GPS-RO observation type (:refractivity, :bending_angle)
- `earth_radius::T`: Earth radius for ray path calculations
- `atmospheric_model::Symbol`: Atmospheric model for ray tracing (:spherical, :ellipsoidal)
- `integration_method::Symbol`: Numerical integration method (:simpson, :trapezoid, :gauss)
"""
struct GPSROOperator{T<:AbstractFloat} <: NonlinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    observation_type::Symbol
    earth_radius::T
    atmospheric_model::Symbol
    integration_method::Symbol
    
    function GPSROOperator{T}(config::ForwardOperatorConfig{T},
                             obs_type::Symbol = :refractivity) where T
        new{T}(
            config, obs_type, T(6371000), :spherical, :simpson
        )
    end
end

"""
    RadarOperator{T<:AbstractFloat}

Forward operator for weather radar observations.

# Fields
- `config::ForwardOperatorConfig{T}`: Operator configuration  
- `radar_type::Symbol`: Radar observation type (:reflectivity, :velocity)
- `radar_wavelength::T`: Radar wavelength in meters
- `beam_width::T`: Radar beam width in degrees
- `hydrometeor_types::Vector{Symbol}`: Hydrometeor types to consider
"""
struct RadarOperator{T<:AbstractFloat} <: NonlinearForwardOperator{T}
    config::ForwardOperatorConfig{T}
    radar_type::Symbol
    radar_wavelength::T
    beam_width::T
    hydrometeor_types::Vector{Symbol}
    
    function RadarOperator{T}(config::ForwardOperatorConfig{T},
                             radar_type::Symbol = :reflectivity) where T
        new{T}(
            config, radar_type, T(0.10), T(1.0),
            [:rain, :snow, :graupel, :hail]
        )
    end
end

"""
    apply_forward_operator(operator::AbstractForwardOperator{T},
                          state::AbstractStateVector{T},
                          observation::AbstractObservation{T}) where T

Apply forward operator to compute observation equivalent value.

# Arguments
- `operator::AbstractForwardOperator{T}`: Forward operator to apply
- `state::AbstractStateVector{T}`: Model state vector  
- `observation::AbstractObservation{T}`: Observation defining location and type

# Returns
- `T`: Computed observation equivalent (H(x))

# Examples
```julia
temp_op = TemperatureOperator{Float64}(config, pressure_levels)
obs_equivalent = apply_forward_operator(temp_op, model_state, temperature_obs)
```
"""
function apply_forward_operator(operator::TemperatureOperator{T},
                               state::AbstractStateVector{T},
                               observation::AbstractObservation{T}) where T
    
    # Get observation location
    obs_location = observation.location
    
    # Compute or retrieve interpolation weights
    location_key = (obs_location.latitude, obs_location.longitude, obs_location.pressure)
    
    if operator.config.cache_interpolation_weights && haskey(operator.cached_weights, location_key)
        weights = operator.cached_weights[location_key]
    else
        weights = compute_interpolation_weights(operator, obs_location, state)
        if operator.config.cache_interpolation_weights
            operator.cached_weights[location_key] = weights
        end
    end
    
    if !weights.is_valid
        return T(NaN)  # Outside valid domain
    end
    
    # Extract temperature field from state vector
    temp_field = get_temperature_field(state)
    
    # Perform interpolation
    interpolated_temp = zero(T)
    for (idx, weight) in zip(weights.grid_indices, weights.weights)
        interpolated_temp += weight * temp_field[idx]
    end
    
    # Apply lapse rate correction if needed
    if operator.apply_lapse_rate_correction
        model_elevation = get_model_elevation_at_location(obs_location, state)
        elevation_diff = obs_location.altitude - model_elevation
        
        # Standard atmospheric lapse rate: 6.5 K/km
        lapse_rate = T(0.0065)  # K/m
        interpolated_temp -= lapse_rate * elevation_diff
    end
    
    return interpolated_temp
end

function apply_forward_operator(operator::HumidityOperator{T},
                               state::AbstractStateVector{T},
                               observation::AbstractObservation{T}) where T
    
    obs_location = observation.location
    location_key = (obs_location.latitude, obs_location.longitude, obs_location.pressure)
    
    # Compute interpolation weights
    if operator.config.cache_interpolation_weights && haskey(operator.cached_weights, location_key)
        weights = operator.cached_weights[location_key]
    else
        weights = compute_interpolation_weights(operator, obs_location, state)
        if operator.config.cache_interpolation_weights
            operator.cached_weights[location_key] = weights
        end
    end
    
    if !weights.is_valid
        return T(NaN)
    end
    
    # Extract humidity and temperature fields
    humidity_field = get_humidity_field(state)
    
    # Interpolate humidity
    interpolated_humidity = zero(T)
    for (idx, weight) in zip(weights.grid_indices, weights.weights)
        interpolated_humidity += weight * humidity_field[idx]
    end
    
    # Convert humidity type if necessary
    if operator.humidity_type == :relative
        # Convert specific humidity to relative humidity
        temp_field = get_temperature_field(state)
        interpolated_temp = zero(T)
        for (idx, weight) in zip(weights.grid_indices, weights.weights)
            interpolated_temp += weight * temp_field[idx]
        end
        
        # Compute saturation vapor pressure and relative humidity
        pressure_pa = obs_location.pressure
        relative_humidity = specific_to_relative_humidity(
            interpolated_humidity, interpolated_temp, pressure_pa
        )
        
        # Apply saturation check
        if operator.saturation_check
            relative_humidity = min(relative_humidity, T(100.0))  # Max 100%
        end
        
        return relative_humidity
    elseif operator.humidity_type == :mixing_ratio
        # Convert specific humidity to mixing ratio: r = q / (1 - q)
        return interpolated_humidity / (1 - interpolated_humidity)
    else
        # Return specific humidity directly
        return interpolated_humidity
    end
end

function apply_forward_operator(operator::WindOperator{T},
                               state::AbstractStateVector{T},
                               observation::AbstractObservation{T}) where T
    
    obs_location = observation.location
    location_key = (obs_location.latitude, obs_location.longitude, obs_location.pressure)
    
    # Compute interpolation weights
    if operator.config.cache_interpolation_weights && haskey(operator.cached_weights, location_key)
        weights = operator.cached_weights[location_key]
    else
        weights = compute_interpolation_weights(operator, obs_location, state)
        if operator.config.cache_interpolation_weights
            operator.cached_weights[location_key] = weights
        end
    end
    
    if !weights.is_valid
        return T(NaN)
    end
    
    # Extract wind fields
    u_field = get_u_wind_field(state)
    v_field = get_v_wind_field(state)
    
    # Interpolate wind components
    interpolated_u = zero(T)
    interpolated_v = zero(T)
    
    for (idx, weight) in zip(weights.grid_indices, weights.weights)
        interpolated_u += weight * u_field[idx]
        interpolated_v += weight * v_field[idx]
    end
    
    # Apply coordinate transformation if needed
    if operator.coordinate_system == :geographic && operator.rotation_angles !== nothing
        # Transform from grid coordinates to geographic coordinates
        lat_idx, lon_idx = get_grid_indices(obs_location)
        rotation_angle = operator.rotation_angles[lat_idx, lon_idx]
        
        cos_rot = cos(rotation_angle)
        sin_rot = sin(rotation_angle)
        
        u_geo = cos_rot * interpolated_u - sin_rot * interpolated_v
        v_geo = sin_rot * interpolated_u + cos_rot * interpolated_v
        
        interpolated_u, interpolated_v = u_geo, v_geo
    end
    
    # Return requested wind component
    if operator.wind_component == :u
        return interpolated_u
    elseif operator.wind_component == :v
        return interpolated_v
    elseif operator.wind_component == :speed
        return sqrt(interpolated_u^2 + interpolated_v^2)
    elseif operator.wind_component == :direction
        # Wind direction in degrees (meteorological convention)
        direction = atan(interpolated_u, interpolated_v) * 180 / π
        return mod(direction + 180, 360)  # Convert to 0-360° range
    else
        throw(ArgumentError("Unknown wind component: $(operator.wind_component)"))
    end
end

function apply_forward_operator(operator::PressureOperator{T},
                               state::AbstractStateVector{T},
                               observation::AbstractObservation{T}) where T
    
    obs_location = observation.location
    
    # Get pressure field from state
    pressure_field = get_pressure_field(state)
    
    if operator.pressure_type == :surface
        # Surface pressure interpolation
        location_key = (obs_location.latitude, obs_location.longitude, T(0))
        
        if operator.config.cache_interpolation_weights && haskey(operator.cached_weights, location_key)
            weights = operator.cached_weights[location_key]
        else
            weights = compute_surface_interpolation_weights(operator, obs_location, state)
            if operator.config.cache_interpolation_weights
                operator.cached_weights[location_key] = weights
            end
        end
        
        if !weights.is_valid
            return T(NaN)
        end
        
        # Interpolate surface pressure
        interpolated_pressure = zero(T)
        for (idx, weight) in zip(weights.grid_indices, weights.weights)
            # Extract surface level (assume last vertical index is surface)
            surface_idx = CartesianIndex(idx[1], idx[2], size(pressure_field, 3))
            interpolated_pressure += weight * pressure_field[surface_idx]
        end
        
        # Apply terrain correction if observation is not at model surface elevation
        if operator.terrain_elevation !== nothing
            model_elevation = interpolate_terrain_elevation(obs_location, operator.terrain_elevation)
            elevation_diff = obs_location.altitude - model_elevation
            
            if abs(elevation_diff) > 10.0  # Only apply correction for significant differences
                # Hydrostatic adjustment: dp/dz = -ρg ≈ -p/(RT) * g
                temp_field = get_temperature_field(state)
                interpolated_temp = zero(T)
                for (idx, weight) in zip(weights.grid_indices, weights.weights)
                    surface_idx = CartesianIndex(idx[1], idx[2], size(temp_field, 3))
                    interpolated_temp += weight * temp_field[surface_idx]
                end
                
                # Apply hydrostatic correction
                R_dry = T(287.05)  # J/(kg·K)
                g = T(9.80665)     # m/s²
                
                pressure_correction = exp(-g * elevation_diff / (R_dry * interpolated_temp))
                interpolated_pressure *= pressure_correction
            end
        end
        
        return interpolated_pressure
        
    elseif operator.pressure_type == :sea_level
        # Sea level pressure reduction
        surface_pressure = apply_forward_operator(
            PressureOperator{T}(operator.config, :surface), state, observation
        )
        
        # Standard sea level pressure reduction
        return reduce_to_sea_level_pressure(surface_pressure, obs_location.altitude)
        
    else
        # Station pressure (same as surface pressure)
        return apply_forward_operator(
            PressureOperator{T}(operator.config, :surface), state, observation
        )
    end
end

function apply_forward_operator(operator::GPSROOperator{T},
                               state::AbstractStateVector{T},
                               observation::AbstractObservation{T}) where T
    
    if operator.observation_type == :refractivity
        return compute_gps_refractivity(operator, state, observation)
    elseif operator.observation_type == :bending_angle  
        return compute_gps_bending_angle(operator, state, observation)
    else
        throw(ArgumentError("Unknown GPS-RO observation type: $(operator.observation_type)"))
    end
end

function apply_forward_operator(operator::RadarOperator{T},
                               state::AbstractStateVector{T}, 
                               observation::AbstractObservation{T}) where T
    
    if operator.radar_type == :reflectivity
        return compute_radar_reflectivity(operator, state, observation)
    elseif operator.radar_type == :velocity
        return compute_radar_velocity(operator, state, observation)
    else
        throw(ArgumentError("Unknown radar observation type: $(operator.radar_type)"))
    end
end

"""
    apply_tangent_linear(operator::AbstractForwardOperator{T},
                        state::AbstractStateVector{T},
                        perturbation::AbstractStateVector{T},
                        observation::AbstractObservation{T}) where T

Apply tangent linear operator: δy = H'(x) δx.

# Arguments
- `operator::AbstractForwardOperator{T}`: Forward operator
- `state::AbstractStateVector{T}`: Reference state x
- `perturbation::AbstractStateVector{T}`: State perturbation δx
- `observation::AbstractObservation{T}`: Observation defining location

# Returns  
- `T`: Tangent linear result δy
"""
function apply_tangent_linear(operator::LinearForwardOperator{T},
                             state::AbstractStateVector{T},
                             perturbation::AbstractStateVector{T},
                             observation::AbstractObservation{T}) where T
    
    # For linear operators, tangent linear is the same as forward operator
    return apply_forward_operator(operator, perturbation, observation)
end

function apply_tangent_linear(operator::HumidityOperator{T},
                             state::AbstractStateVector{T}, 
                             perturbation::AbstractStateVector{T},
                             observation::AbstractObservation{T}) where T
    
    # Tangent linear for humidity operator (considering temperature dependency)
    obs_location = observation.location
    
    # Get interpolation weights
    location_key = (obs_location.latitude, obs_location.longitude, obs_location.pressure)
    if haskey(operator.cached_weights, location_key)
        weights = operator.cached_weights[location_key]
    else
        weights = compute_interpolation_weights(operator, obs_location, state)
    end
    
    if !weights.is_valid
        return T(NaN)
    end
    
    # For specific humidity, tangent linear is linear
    if operator.humidity_type == :specific
        humidity_pert_field = get_humidity_field(perturbation)
        
        result = zero(T)
        for (idx, weight) in zip(weights.grid_indices, weights.weights)
            result += weight * humidity_pert_field[idx]
        end
        
        return result
        
    elseif operator.humidity_type == :relative
        # Tangent linear for relative humidity includes temperature dependency
        # ∂RH/∂q and ∂RH/∂T terms
        
        # This is simplified - full implementation would need proper derivatives
        humidity_pert_field = get_humidity_field(perturbation)
        
        result = zero(T)  
        for (idx, weight) in zip(weights.grid_indices, weights.weights)
            result += weight * humidity_pert_field[idx]
        end
        
        # Scale by approximate conversion factor
        return result * T(100.0)  # Simplified conversion
    else
        # Mixing ratio tangent linear
        humidity_pert_field = get_humidity_field(perturbation)
        humidity_field = get_humidity_field(state)
        
        result = zero(T)
        for (idx, weight) in zip(weights.grid_indices, weights.weights)
            q = humidity_field[idx]
            dq = humidity_pert_field[idx] 
            # d(q/(1-q))/dq = 1/(1-q)²
            result += weight * dq / ((1 - q)^2)
        end
        
        return result
    end
end

"""
    apply_adjoint(operator::AbstractForwardOperator{T},
                 state::AbstractStateVector{T}, 
                 obs_perturbation::T,
                 observation::AbstractObservation{T}) where T -> AbstractStateVector{T}

Apply adjoint operator: δx* = H'ᵀ(x) δy*.

# Arguments  
- `operator::AbstractForwardOperator{T}`: Forward operator
- `state::AbstractStateVector{T}`: Reference state x
- `obs_perturbation::T`: Observation space perturbation δy*
- `observation::AbstractObservation{T}`: Observation defining location

# Returns
- `AbstractStateVector{T}`: Adjoint result δx*
"""
function apply_adjoint(operator::TemperatureOperator{T},
                      state::AbstractStateVector{T},
                      obs_perturbation::T,
                      observation::AbstractObservation{T}) where T
    
    # Create result state vector
    result_state = StateVector(state.config)
    
    obs_location = observation.location
    location_key = (obs_location.latitude, obs_location.longitude, obs_location.pressure)
    
    # Get interpolation weights
    if haskey(operator.cached_weights, location_key)
        weights = operator.cached_weights[location_key]
    else
        weights = compute_interpolation_weights(operator, obs_location, state)
    end
    
    if !weights.is_valid
        return result_state  # Return zero state
    end
    
    # Apply adjoint (transpose of interpolation)
    temp_field = get_temperature_field(result_state)
    
    for (idx, weight) in zip(weights.grid_indices, weights.weights)
        temp_field[idx] += weight * obs_perturbation
    end
    
    return result_state
end

# Helper functions for field extraction and interpolation

"""
    compute_interpolation_weights(operator::AbstractForwardOperator{T},
                                location::ObservationLocation{T},
                                state::AbstractStateVector{T}) where T

Compute interpolation weights for a given observation location.
"""
function compute_interpolation_weights(operator::AbstractForwardOperator{T},
                                     location::ObservationLocation{T},
                                     state::AbstractStateVector{T}) where T
    
    grid_config = state.config  # Assume state has grid configuration
    
    # Convert geographic coordinates to grid indices
    grid_coords = geographic_to_grid_coordinates(
        location.latitude, location.longitude, grid_config
    )
    
    # Handle vertical coordinate
    if operator.config.vertical_coordinate == :pressure
        vertical_coord = location.pressure
    else
        vertical_coord = location.altitude
    end
    
    # Compute interpolation weights based on method
    if operator.config.interpolation_method == BILINEAR_INTERPOLATION
        return compute_bilinear_weights(grid_coords, vertical_coord, state)
    elseif operator.config.interpolation_method == BICUBIC_INTERPOLATION
        return compute_bicubic_weights(grid_coords, vertical_coord, state)
    else
        # Default to bilinear
        return compute_bilinear_weights(grid_coords, vertical_coord, state)
    end
end

function compute_bilinear_weights(grid_coords::Tuple{T,T}, vertical_coord::T,
                                 state::AbstractStateVector{T}) where T
    
    i_float, j_float = grid_coords
    grid_size = state.config.grid_size
    nx, ny, nz = grid_size
    
    # Check bounds
    if i_float < 1 || i_float > nx || j_float < 1 || j_float > ny
        return InterpolationWeights{T}(
            CartesianIndex{3}[], T[], false
        )
    end
    
    # Grid indices and weights for horizontal interpolation
    i1, i2 = floor(Int, i_float), ceil(Int, i_float)
    j1, j2 = floor(Int, j_float), ceil(Int, j_float)
    
    wi = i_float - i1
    wj = j_float - j1
    
    # Find vertical level (simplified pressure coordinate)
    pressure_levels = get_pressure_levels(state)
    k1 = find_pressure_level(vertical_coord, pressure_levels)
    k2 = min(k1 + 1, nz)
    
    if k1 == k2
        wk = T(0)
    else
        wk = (vertical_coord - pressure_levels[k1]) / 
             (pressure_levels[k2] - pressure_levels[k1])
    end
    
    # Bilinear interpolation weights (8 points for 3D)
    indices = CartesianIndex{3}[]
    weights = T[]
    
    for k in [k1, k2], j in [j1, j2], i in [i1, i2]
        # Ensure indices are within bounds
        i_safe = clamp(i, 1, nx)
        j_safe = clamp(j, 1, ny)  
        k_safe = clamp(k, 1, nz)
        
        push!(indices, CartesianIndex(i_safe, j_safe, k_safe))
        
        # Compute trilinear weight
        wi_local = (i == i1) ? (1 - wi) : wi
        wj_local = (j == j1) ? (1 - wj) : wj  
        wk_local = (k == k1) ? (1 - wk) : wk
        
        weight = wi_local * wj_local * wk_local
        push!(weights, weight)
    end
    
    return InterpolationWeights{T}(indices, weights, true)
end

function compute_bicubic_weights(grid_coords::Tuple{T,T}, vertical_coord::T,
                                state::AbstractStateVector{T}) where T
    # Simplified bicubic - full implementation would use 4x4x4 stencil
    return compute_bilinear_weights(grid_coords, vertical_coord, state)
end

# Field extraction functions (would be implemented based on StateVector structure)

function get_temperature_field(state::AbstractStateVector{T}) where T
    # Extract temperature field from state vector
    # This would depend on the specific StateVector implementation
    return reshape(state.variables.t, state.config.grid_size)
end

function get_humidity_field(state::AbstractStateVector{T}) where T
    return reshape(state.variables.q, state.config.grid_size)
end

function get_u_wind_field(state::AbstractStateVector{T}) where T
    return reshape(state.variables.u, state.config.grid_size)
end

function get_v_wind_field(state::AbstractStateVector{T}) where T
    return reshape(state.variables.v, state.config.grid_size)
end

function get_pressure_field(state::AbstractStateVector{T}) where T
    return reshape(state.variables.prse, state.config.grid_size)
end

function get_pressure_levels(state::AbstractStateVector{T}) where T
    # Return model pressure levels
    return collect(range(T(100000), T(1000), length=state.config.grid_size[3]))
end

# Physical conversion functions

function specific_to_relative_humidity(specific_humidity::T, temperature::T, pressure::T) where T
    # Convert specific humidity to relative humidity
    # RH = 100 * (q * pressure) / (0.622 * e_sat(T))
    
    # Saturation vapor pressure (simplified Magnus formula)
    e_sat = T(611.2) * exp(17.67 * (temperature - 273.15) / (temperature - 29.65))
    
    # Mixing ratio from specific humidity
    mixing_ratio = specific_humidity / (1 - specific_humidity)
    
    # Relative humidity
    return 100 * mixing_ratio * pressure / (T(0.622) * e_sat)
end

function reduce_to_sea_level_pressure(surface_pressure::T, elevation::T) where T
    # Standard atmosphere sea level pressure reduction
    # Using barometric formula: p_0 = p * exp(gz/RT)
    
    g = T(9.80665)        # m/s²
    R = T(287.05)         # J/(kg·K)  
    T_std = T(288.15)     # Standard temperature (K)
    
    return surface_pressure * exp(g * elevation / (R * T_std))
end

# GPS-RO forward operators (simplified implementations)

function compute_gps_refractivity(operator::GPSROOperator{T}, 
                                 state::AbstractStateVector{T},
                                 observation::AbstractObservation{T}) where T
    # Simplified GPS refractivity computation
    # Full implementation would integrate along ray path
    return T(300.0)  # Placeholder N-units
end

function compute_gps_bending_angle(operator::GPSROOperator{T},
                                  state::AbstractStateVector{T}, 
                                  observation::AbstractObservation{T}) where T
    # Simplified GPS bending angle computation  
    return T(0.001)  # Placeholder radians
end

# Radar forward operators (simplified implementations)

function compute_radar_reflectivity(operator::RadarOperator{T},
                                   state::AbstractStateVector{T},
                                   observation::AbstractObservation{T}) where T
    # Simplified radar reflectivity computation
    return T(35.0)  # Placeholder dBZ
end

function compute_radar_velocity(operator::RadarOperator{T}, 
                               state::AbstractStateVector{T},
                               observation::AbstractObservation{T}) where T
    # Simplified radar velocity computation
    return T(10.0)  # Placeholder m/s
end

# Helper functions

function find_pressure_level(target_pressure::T, pressure_levels::Vector{T}) where T
    # Find pressure level index for interpolation
    for k in 1:length(pressure_levels)-1
        if pressure_levels[k] >= target_pressure >= pressure_levels[k+1]
            return k
        end
    end
    return length(pressure_levels)
end

function get_model_elevation_at_location(location::ObservationLocation{T}, 
                                        state::AbstractStateVector{T}) where T
    # Get model terrain elevation at observation location
    return T(0.0)  # Placeholder - would interpolate from terrain field
end

function interpolate_terrain_elevation(location::ObservationLocation{T}, 
                                      terrain::Array{T,2}) where T
    # Interpolate terrain elevation to observation location
    return T(0.0)  # Placeholder
end

function get_grid_indices(location::ObservationLocation{T}) where T
    # Convert geographic location to grid indices
    return (1, 1)  # Placeholder
end

function compute_surface_interpolation_weights(operator::PressureOperator{T},
                                             location::ObservationLocation{T},
                                             state::AbstractStateVector{T}) where T
    # Compute weights for surface-level interpolation
    grid_coords = (T(1), T(1))  # Placeholder
    return compute_bilinear_weights(grid_coords, T(0), state)
end

"""
    create_operator_suite(config::ForwardOperatorConfig{T}) where T

Create a complete suite of forward operators for all observation types.

# Arguments
- `config::ForwardOperatorConfig{T}`: Configuration for all operators

# Returns
- `Dict{Symbol,AbstractForwardOperator{T}}`: Dictionary of operators keyed by observation type

# Examples
```julia
config = ForwardOperatorConfig{Float64}()
operators = create_operator_suite(config)

# Use temperature operator
temp_equivalent = apply_forward_operator(
    operators[:temperature], model_state, temperature_obs
)
```
"""
function create_operator_suite(config::ForwardOperatorConfig{T}) where T
    pressure_levels = collect(range(T(100000), T(1000), length=64))
    
    operators = Dict{Symbol,AbstractForwardOperator{T}}(
        :temperature => TemperatureOperator{T}(config, pressure_levels),
        :humidity => HumidityOperator{T}(config, :specific),
        :relative_humidity => HumidityOperator{T}(config, :relative),
        :wind_u => WindOperator{T}(config, :u, :geographic),
        :wind_v => WindOperator{T}(config, :v, :geographic),
        :wind_speed => WindOperator{T}(config, :speed, :geographic),
        :surface_pressure => PressureOperator{T}(config, :surface),
        :sea_level_pressure => PressureOperator{T}(config, :sea_level),
        :gps_refractivity => GPSROOperator{T}(config, :refractivity),
        :gps_bending_angle => GPSROOperator{T}(config, :bending_angle),
        :radar_reflectivity => RadarOperator{T}(config, :reflectivity),
        :radar_velocity => RadarOperator{T}(config, :velocity)
    )
    
    return operators
end

end # module ForwardOperators